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We present Monte Carlo simulations of the Sivers effect in the polarized Drell-Yan pp T —> p + p~X 
process at the center-of-mass energy y/s = 200 GeV reachable at the Relativistic Heavy-Ion Collider 
(RHIC) of BNL. We use two different parametrizations for the Sivers function, one deduced from 
the analysis of Semi-Inclusive Deep-Inelastic Scattering (SIDIS) data at much lower energies, and 
another one constrained by the RHIC data for the pp^ —> tvX process at the same energy. For 
a given luminosity of 10 32 cm -2 sec -1 , we explore the necessary conditions to reach a statistical 
accuracy that allows to extract unambiguous information on the structure of the Sivers function. In 
particular, we consider the feasibility of the test on its predicted universality property of changing 
sign when switching from SIDIS to Drell-Yan processes. 

PACS numbers: 13.88+e,13.85.-t,13.85.Qk 


I. INTRODUCTION 

Azimuthal asymmetries in hard collisions involving (polarized) hadrons represent a formidable testground for Quan¬ 
tum ChromoDynamics (QCD) in the nonperturbative regime. Since almost thirty years, data have been collected 
for hadron-hadron collisions, and recently also for semi-inclusive 7 *-hadron processes in the regime of Deep-Inelastic 
Scattering (DIS). Large asymmetries were observed in the azimuthal distribution of final-state products (with respect 
to the normal to the production plane), particularly when flipping the transverse polarization of one hadron involved 
in the initial or final state: the socalled Single-Spin Asymmetries (SSA). Examples of such SSA are found for the 
pp —> A^X process [1] at forward rapidity, for the pp 1 ' —> tyX process [2-5] again at forward rapidity, and for the semi- 
inclusive hadron production 7 *A) —> 7 tX, where A is the proton [6-9] or the deuteron [10]. Apart from the last one, in 
all other cases SSA up to 40% were detected which were totally unexpected, since they cannot be easily accommodated 
in a consistent manner within the perturbative QCD in the collinear massless approximation [11]; moreover, they seem 
to persist also at higher energies typical of the collider regime [3-5], which also contradicts QCD expectations. The 
same conclusion holds also for unpolarized Drell-Yan experiments at high energy like 7 tA —> p + p~ X [12-15], with 
A=p,d, W, and pp —» p + p~ X [15], where the violation of the Lam-Tung sum rule strongly supports the conjecture 
to go beyond the collinear approximation [14]. 

All this amount of puzzling measurements have triggered an intense theoretical activity, particularly about the idea 
that intrinsic transverse momenta of partons, together with transverse spin degrees of freedom, could be responsible 
for the observed asymmetries. Transverse-Momentum Dependent (TMD) parton distributions and fragmentation 
functions have been introduced and linked to measurable asymmetries in the leading-twist cross sections of Semi- 
Inclusive DIS (SIDIS), Drell-Yan process, semi-inclusive hadron-lradron collision, and e + e - annihilation [16-25]. As 
for parton distributions, the prototype of such TMD functions is the Sivers function [20], which has a probabilistic 
interpretation: it describes how the distribution of unpolarized quarks is distorted by the transverse polarization 
of the parent hadron. Using the notations recommended in Ref. [26], the Sivers function f^ T can be extracted by 
measuring the socalled Sivers effect in hadron-hadron collisions or SIDIS processes, i.e. an asymmetric distribution 
of the final-state products in the azimuthal angle defined by the mixed product p T x P ■ S T , where P is the nucleon 
momentum and p T , S T are the transverse components of the parton momentum and of the nucleon spin with respect 
to the direction of P in the infinite momentum frame. Time-reversal invariance would forbid such correlation if 
there were no initial/final-state interactions in the considered collision/SIDIS process, respectively. Therefore, f^ T is 
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conventionally named a ” naive time-reversal-odd” distribution. The interactions must imply an interference between 
different helicity states of the target nucleon [27, 28]; consequently, the correlation between p T and S T is possible 
only for a nonvanishing orbital angular momentum of the partons. Then, extraction of f^ T from data allows to study 
the orbital motion of hidden confined partons; better, it contains information on their spatial distribution [29], and 
it offers a natural link between microscopic properties of confined elementary constituents and hadronic measurable 
quantities, such as the nucleon anomalous magnetic moment [30]. 

In QCD, the necessary interactions can be naturally identified with the multiple exchange of soft gluons contained in 
the gauge link operator, which grants the color gauge-invariant definition of TMD distributions [31, 32]. However, the 
whole picture relies on the proof of a suitable factorization theorem at small trasverse momenta for the process at hand. 
At present, QCD factorization proofs have been established for e + e - annihilations [17], for Drell-Yan processes [19], 
and, more recently, for SIDIS processes [33, 34], including also naive T-odd contributions. The related universality 
of TMD functions has been carefully discussed in Refs. [31, 33, 35]. It turns out that the Sivers function displays 
the very interesting property of changing sign when going from the SIDIS to the Drell-Yan process, due to a peculiar 
feature of its gauge link operator under the time-reversal operation [31]. This interesting prediction has stimulated 
intense experimental and phenomenological activities to link the Sivers effect recently measured at HERMES [7, 8] 
with processes happening at RHIC, where data are being taken for polarized pp collisions [4]. In particular, three 
different parametrizations of f^p [36-38] have been extracted from the HERMES data (and found compatible also 
with the recent COMPASS data [10]), and have been used then to make predictions for SSA at RHIC (see also the 
more recent analysis of Ref. [39]; for a comparison among the various approaches, see also Ref. [40]). 

In view of the foreseen upgrade of RHIC detector and luminosity (RHIC II), we will consider here the specific Drell- 
Yan process pp^ —> p + X. The leading-twist polarized part of the cross section contains two terms that produce 
interesting SSA with azimuthal distinct behaviours [41]. In a previous paper [42], we analyzed the term weighted by 
sin {(j) + (f>s), with (j> and <j>s the azimuthal orientations of the final lepton plane and of the proton polarization with 
respect to the reaction plane; it leads to the extraction of another interesting naive T-odd TMD distribution, the 
Boer-Mulders h^, which is most likely responsible for the above mentioned violation of the Lam-Tung sum rule [41]. 
In Ref. [42], we considered the Drell-Yan process pp^ —> fi + p~X at the kinematics of interest for the High Energy 
Storage Ring (HESR) project at GSI [43, 44] and we numerically simulated the SSA with a Monte Carlo in order to 
explore the minimal conditions required for an unambiguous extraction of hj~. Here, we follow the same approach to 
isolate the other term weighted by sin((f> — <j>s), which contains the convolution of fyp with the standard unpolarized 
parton distribution f\. The Monte Carlo will be applied to the pp"' —> p + X process at yfs = 200 GeV at the 
RHIC II luminosity (at least 10 32 cm -2 sec -1 ). The SSA will be numerically simulated using as input both the 
parametrization of Ref. [36] and a new high-energy parametrization of f^ T constrained by recent RHIC data on SSA 
for the pp'' —> 7 tX process at the same yfs = 200 GeV [4]. The goal is to explore the sensitivity of the simulated 
asymmetry to different input parametrizations, as well as to directly verify, within the reached statistical accuracy, 
the predicted sign change of the Sivers function between SIDIS and Drell-Yan [31]. 

In Sec. II, the general formalism and details of the numerical simulation are briefly reviewed. In Sec. Ill, we discuss 
the input parametrizations. In Sec. IV, results are presented. Finally, in Sec. V some conclusions are drawn. 



II. THEORETICAL FRAMEWORK AND NUMERICAL SIMULATIONS 

In a Drell-Yan process, a lepton with momentum k\ and an antilepton with momentum &2 (with k 2 , 2 ) ~ 0) are 
produced from the collision of two hadrons with momentum Pi, mass M\, spin Si, and P 2 , M 2 , S 2 , respectively (with 
^ 1 ( 2 ) = -^ 1 ( 2 )> ^ 1 ( 2 ) = — 1) Pi( 2 ) • <S 1 ( 2 ) = 0). The center-of-mass (cm) square energy available is s = (Pi + P 2) 2 
and the invariant mass of the final lepton pair is given by the time-like momentum transfer q 2 = M 2 = ( k\ + k 2 ) 2 - 
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In the kinematical regime where M 2 ,s —■> oo, while keeping the ratio 0 < r = M 2 /s < 1 limited, the lepton pair 
can be assumed to be produced from the elementary annihilation of a parton and an antiparton with momenta pi 
and P 2 , respectively. If P+ and P 2 _ are the dominant light-cone components of hadron momenta in this regime, 
then the partons are approximately collinear with the parent hadrons and carry the light-cone momentum fractions 
0 < x\ = Pi/P^, X 2 = P 2 IP -2 — 1) with q + = Pi , q~ = pi by momentum conservation [41]. As already 
anticipated in Sec. I, a key issue is the extraction of TMD parton distributions; this requires the cross section to 
be kept differential in the transverse momentum of the final lepton pair, q T , which is bounded by the momentum 
conservation q T = p\ T + p 2 T to each intrinsic transverse components pi T of the parton momentum pi with respect to 
the direction defined by the corresponding hadron momentum Pi. If q T ^ 0 the annihilation direction is not known. 
Hence, it is convenient to select the socalled Collins-Soper frame [45] (see fig. 1), where 
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The final lepton pair is detected in the solid angle (0, </>), where 6 is defined in Fig. 1 and <j> (and all other azimuthal 
angles) is measured in a plane perpendicular to z, t, but containing h. 

The full expression of the leading-twist differential cross section for the ~ > l + l~X process can be written 

as [41] 
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where a is the fine structure constant, dQ, = sin 6d8d<j>, e/ is the charge of the parton with flavor /, <pSi is the azimuthal 
angle of the transverse spin of hadron i. and 

A (y) = - y + y 2 ^j \ (! + cos 2 e ) B (y ) = vi 1 -y) C = \ sin 2 9 . (3) 

The TMD functions f{(H), hi^(H), describe the distributions of unpolarized and transversely polarized partons in 
an unpolarized hadron H, respectively, while .fy/(H), h^jf (H), have a similar interpretation but for transversely 
polarized hadrons . The transversity h{ describes transversely polarized partons in transversely polarized hadrons. 
Each one of these distributions for a parton / is convoluted with its antiparton partner / according to 


T [Af (Hi) A*(H 2 )\ 
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In previous papers, we analyzed the SSA generated by the azimuthal dependences cos 20 and sin(0 + <ps 2 ) in 
Eq. (2) [42], as well as the double-polarized Drell-Yan process [46]. A combined measurement of these SSA allows 
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to completely determine the unknown transversity hi and Boer-Mulders function , which could be responsible for 
the well known violation of the Lam-Tung sum rule [41] in unpolarized Drell-Yan data [12- 14] (see also Ref. [47], 
and references therein, for a recent discussion on a parallel with QCD vacuum effects). We set up a Monte Carlo 
simulation of Drell-Yan processes involving unpolarized antiproton beams and transversely polarized protons for 
various kinematic scenarios at the HESR at GSI, namely for 30 < s < 200 GeV 2 with an antiproton beam energy 
of 15 GeV and the socalled asymmetric collider mode (proton beams of 3 GeV) or fixed target mode (proton fixed 
targets). Special focus was put on the range 4 < M < 9 GeV for the lepton invariant mass, since it does not overlap 
with the charmonium and bottonium resonances (where the elementary annihilation does not necessarily proceed 
through a simple intermediate virtual photon) and higher-twist corrections should be suppressed justifying the simple 
approach of Eq. (2) based on the parton model [42]. Here, we will concentrate on the term weighted by sin(<(> — <^g 2 ) 
in Eq. (2) for Hi = p and Hi = p\ and we will consider the related SSA at \fs = 200 GeV reachable at RHIC [48]. 
Most of the technical details of the simulation are mutuated from our previous works; hence, we will heavily refer to 
Refs. [42, 46] in the following. 

The Monte Carlo events have been generated by the following cross section [42]: 


da 

dfldxidx2dq T 


K - A(q T , xi,X 2 , M) F{xi,x 2 ) Y, d(q T ,x i,x 2 ) S*(0,0, <j> S2 ) ■ 
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It means that in Eq. (2) we assume a factorized transverse-momentum dependence in each parton distribution such 
as to break the convolution F, leading to the product A F. The function A is parametrized as [42] 
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where a(x F ,M), b(x F ,M), are parametric polynomials given in Appendix A of Ref. [14] with x F = Xi — X 2 and 
q T = \q T \ (see also the more recent Ref. [49]). It is normalized as 


J dq T A(q T ,xi,x 2 ,M) = 1 . 


(7) 


Actually, the Drell-Yan events studied in Ref. [14] were produced for 7r — p collisions; however, the same analysis, 
repeated for p — p and p — p collisions [15], gives a similar distribution for q T not very close to 0 and not much larger 
than 3 GeV/c. Here, we will adopt two different cuts on q T depending on the input parametrization for the Sivers 
function /-(y (see next Sec. Ill), namely 1 < q T < 3 GeV/c and 0.1 < q T <2 GeV/c. Anyway, the average transverse 
momentum turns out to be (q T ) > 1 GeV/c, i.e. much bigger than parton intrinsic transverse momenta induced by 
confinement. 

The latter observation implies that sizeable QCD corrections affect the simple parton model picture of Eq. (2). Their 
influence on the q T distribution is effectively contained in the phenomenological parametrization of Eq. (6). However, 
there are also other well known corrections [50] coming from the resummation of leading logarithms at any order in 
the strong coupling constant a s , and from the inclusion of diagrams at first order in a s involving the qq fusion or the 
gy Compton mechanisms. The first group, usually named leading-log approximation (LLA), introduces a logarithmic 
dependence on the scale M 2 inside the various parameters entering the parton distributions [51] contained in Eq. (5), 
such that it would determine their DGLAP evolution. However, the range of M values here explored is close to the 
one of Refs. [14, 15], where the parametrization of A,F, and Cj in Eq. (5) was deduced assuming M-independent 
parton distributions. Moreover, as it will be shown in the next Sec. IV, most of the events concentrate around the 
average (x) ~ 0.1, where the effects of evolution are almost vanishing. Therefore, similarly to Ref. [42, 46] we take 

2 

F(xI,x 2 ) = Y e f fl(xi'J/Hi) f{(x 2 ;f/H 2 ) + (f <-> f) , (8) 

which represents the azimuthally symmetric unpolarized part of Eq. (2) that has been factorized out for convenience. 
As previously stressed, the unpolarized distribution f{ (x) for various flavors / = it, d, s, is parametrized as in Ref. [15]. 

The second group of QCD corrections is named next-to-leading-log approximation (NLLA), and it is responsible 
for the well known K factor in Eq. (5). The K factor is roughly independent on x F and M 2 but it can grow with 
\/t\ it also depends on the chosen normalization of the parton distributions (for a detailed analysis in this context see 
Ref. [15]). It is a large correction, typically a multiplicative factor in the range 1.5-7 2.5. Here, we will conventionally 
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assume the same value 2.5 adopted in our previous simulations. But we stress that in an azimuthal asymmetry the 
corrections to the cross sections in the numerator and in the denominator should compensate each other. This is 
certainly true for each elementary contribution to the amplitude for SSA, but it is much less obvious for the ratio of 
full differential cross sections. Indeed, the smooth dependence of the SSA on NLLA corrections has been confirmed 
for fully polarized Drell-Yan processes at RHIC cm square energies [52]. 

The whole solid angle ( 9 , cf >) of the final lepton pair in the Collins-Soper frame is randomly distributed in each 
variable. The explicit form for sorting the angular distribution in the Monte-Carlo is [42, 46] 

y; Cj(q T ,x 1,3:2) Si(fl,(ft,<fts 2 ) = 1 + cos 2 9 + u ^ Xl ' x ^ 2 ' qT ^ sin 2 flcos 2 <ft 

i=1 

+ \S 2T \c 4 {q T ,Xi,X 2 ) SliO, </>, (f>s 2 ) ■ (9) 

If quarks were massless, the virtual photon would be only transversely polarized and the angular dependence would 
be described by the functions C\ = S\ — 1 and C 2 = 1, S 2 = cos 2 9. Violations of such azimuthal symmetry induced 
by the function 03 = ^ are due to the longitudinal polarization of the virtual photon and to the fact that quarks 
have an intrinsic transverse momentum distribution, leading to the explicit violation of the socalled Lam-Tung sum 
rule [14], QCD corrections influence v, which in principle depends also on M 2 (see App. A of Ref. [14]). Azimuthal 
cos 2 <j> asymmetries induced by v were simulated in Ref. [42] using the simple parametrization of Ref. [41] and testing 
it against the previous measurement of Ref. [14]. 

The last term in Eq. (9) corresponds to the polarized part of the cross section (2). Since we want to single out just 
the Sivers contribution, we assume that 

Sa{ 9, 4>, (f>s 2 ) = (1 + cos 2 9) sin (0 - (j>s 2 ) ■ ( 10 ) 

Recalling that in Eq. (5) the azimuthally symmetric unpolarized part A(q T , 24 , x 2 , M) F(x\, # 2 ) of the cross section 
has been factorized out, the corresponding coefficient C 4 in Eq. (9) in principle reads 
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where the complete dependence of the involved TMD parton distributions has been made explicit. In the next 
Sec. Ill, we will discuss two different parametrizations for the x and p T dependence of these distributions which allow 
to calculate the convolutions and determine C4. 

Following Refs. [42, 46], the SSA corresponding to the Sivers effect is constructed by dividing the event sample in 
two groups, one for positive values of sin (</> — <j>s 2 ) (U) and another one for negative values ( D ), and taking the ratio 
(U — D)/(U + D ). Data are accumulated only in the x 2 bin, i.e. they are summed over in aq, 9 , and in q T with the 
discussed cutoffs. Contrary to Refs. [42, 46], no cutoff has been applied to the 9 distribution because S 4 in Eq. (10) 
contains the term 1 + cos 2 9. Statistical errors for (U — D)/(U + D) are obtained by making 10 independent repetitions 
of the simulation for each individual case, and then calculating for each X 2 bin the average asymmetry value and the 
variance. We checked that 10 repetitions are a reasonable threshold to have stable numbers, since the results do not 
change significantly when increasing the number of repetitions beyond 6. In a real experiment, the SSA would be 
extracted by taking the ratio between proper differences and sums of cross sections for the four possible combinations 
with the azimuthal angles ±</>, ±</>s 2 , in order to reduce systematic errors. In the Monte Carlo simulation, for each 
(j>g 2 we can simply build the SSA in the </> angle. In an ideal experiment, the two situations would be equivalent. It 
is worth noting that while </>g 2 is fixed in the lab frame, in the Collins-Soper frame of Fig. 1 it is variable, since the h 
axis is directed along q T /g T ; hence, a random distribution in <j>s 2 must be initially extracted in the Monte Carlo. 


III. PARAMETRIZATIONS OF THE SIVERS FUNCTION 

In our previous papers [42, 46], the strategy of the numerical simulation was based on making guesses for the input 
x and p T dependence of the parton distributions, and on trying to determine the minimum number of events required 
to discriminate various SSA produced by very different input guesses. In fact, this would be equivalent to state that 
in this case some analytic information on the structure of these TMD parton distributions could be extracted from 
the SSA measurement. 
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As for the Sivers effect, the situation is different because recently the HERMES collaboration has released new 
SSA data for the SIDIS process on transversely polarized protons [8], which substantially increase the precision of the 
previous data set [7]. As a consequence, three different parametrizations of f^p [36-38] have been extracted from this 
data set and found compatible also with the recent COMPASS data [10]. Moreover, a recent preprint appeared which 
usefully illustrates the differences among the various approaches [40]. At the same time, new data have been collected 
at RHIC [4, 5] on SSA in the pp^ —> irX process, that confirm the observation of large asymmetries at forward rapidity 
of the pion also at the high-energy collider regime (y^ = 200 GeV). Despite this class of hadron-hadron collisions is 
power suppressed and factorization was established in terms of higher-twist correlation operators [53], still the SSA 
receives a contribution from the leading-twist convolution /i ® (g> D\ , where D\ describes the fragmentation of an 

unpolarized quark into the detected 7r. Therefore, analogously to the analysis of Ref. [54] at lower energy [2], we believe 
that the measured SSA in pp^ —> nX processes at \/s = 200 GeV should indirectly constrain the parametrization of 
and, consequently, the ’’strength” of the Sivers effect when these data are ideally interpreted as completely driven 
by the above convolution. 

In our Monte Carlo simulation, we consider two different parametrizations for the one elaborated in Ref. [36] 
and based on the new HERMES [8] and COMPASS [10] data, where the p T dependence is driven by the (p 2 ) extracted 
in a model dependent way from the azimuthal asymmetry of the unpolarized SIDIS cross section (Cahn effect); a new 
high-energy parametrization inspired to the one of Ref. [37] but with a specific p T dependence constrained by the 
pp T TrV'data at s/s = 200 GeV. 


1. The parametrization of Ref. [36] 


Keeping in mind the commonly adopted conventions [26], the expression used in Ref. [36] is 

AV(*.Pt) = ^ N ff/pf(x,PT) 

-2 N f (af + b f)“ f+bf x a, (1 _ x )bf f{( X ,p T ) 

} at ’ bf v ’ pl + M 2 y J 


ay b b / 
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IT {pi) 


(«/ + &f) a/+b/ x a f _ x p e -p 2 T/(p 2 T ) //( x ) 

a a /b b / Pl + M 2 JlK 


( 12 ) 


where M 2 is the mass of the polarized proton, p T = \p T \, and (p 2 ) = 0.25 (GeV/c) 2 is deduced by assuming a Gaussian 
form for the p T dependence of fi in order to reproduce the azimuthal angular dependence of the SIDIS unpolarized 
cross section (Cahn effect). The parameters Nf, af, bf, with f = u,d are extracted by fitting the recent HERMES [8] 
and COMPASS [10] data with a final % 2 per degree of freedom of 1.06 (the negligible contribution of antiquarks in 
the minimization procedure has been traded off for a better precision). They are listed in Tab. I. 


TABLE I: Parameters for the Sivers distribution from Ref. [36] 


quark up 


quark down 


N u 

0.32 ±0.11 

N d 

-1.0 ±0.12 

O'U 

0.29 ±0.35 

a d 

1.16 ±0.47 

bu 

0.53 ±3.58 

b d 

3.77 ±2.59 

A/ 0 2 

0.32 ± 0.25 (GeV/c) 2 




Following the prediction about a sign change of when going from a SIDIS process (as for the HERMES 
analysis) to a Drell-Yan process (as it is considered here), we insert the opposite of Eq. (12) inside Ecp (11), including 
the Gaussian parametrization for f\(x,p T ). The integrals upon the transverse momenta can be evaluated following 
the steps in Sec.VI of Ref. [41]. The net result is 

E/ e 2 Nf i ' af + a b /\ - fi(x 1 ;f/p)xy (1 - x 2 ) bf /l)^;//p T ) 

r | JMoqr _ a/ V _ 

C4 - '* 2r| , 2 + 4 MS E, ^ fiM/p) f!(xr,f/p) + (1 - 2) ’ 

We further simplify this expression by replacing the flavor-dependent product of parton distributions with an average 
product ( fi(xi )) (fi(x 2 )) both in the numerator and in the denominator, in order to reduce the statistical noise related 
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to the parametrization of fi(x). The final expression of C 4 for pp^ collisions, that numerically simulates the Sivers 
effect in our Monte Carlo, becomes 1 


c 4 


4Mo q T 1 
q 2 T + 4Mq 9 
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k + b u r^_ ^ (1 _ ^ + ^ ^+ a ^ d+&d (i _ 


70 

a u o u 
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(14) 


As it is evident from previous formulae, the parametrization ( 12 ) is put in a very convenient form that easily fits the 
asymmetry term c 4 in our Monte Carlo. However, the sometimes poor resolution in determining the parameters forced 
us to select only the central values in Tab. I in order to produce meaningful numerical simulations. The sensitivity 
of the parameters to the HERMES results for the Sivers effect reflects in a more important relative weight of the d 
quark over the u one in the valence x range, with opposite signs for the corresponding normalization Nf, f = u,d. 
As it will be shown in the next Sec. IV, this has two main consequences on the simulation: small SSA are obtained 
for the pp^ — > p + p~X process in the valence x range, where the dd annihilation occurs less frequently than the 
uu one; a significant minimum number of events is necessary in the sample to reduce the statistical error bars and 
make the asymmetry not compatible with zero. From Tab. I, the a u parameter is much smaller than 1. This means 
that in Eq. (14) the u-quark term dominates at small X 2 , leading to a persistence of the Sivers effect even below the 
valence x range. This feature is potentially very relevant at RHIC kinematics, where (x) ~ 0.01. Therefore, for this 
parametrization we have produced also a specific simulation a small X 2 with a finer binning A 2 * 2 , as it will be shown 
in the next Sec. IV. The flavor-independent Lorentzian shape in the p T dependence of Eq. (12) produces a maximum 
asymmetry for q T ~ 1 GeV/c and a rapid decrease for larger values. Consequently, transverse momenta are selected 
in the range 0.1 < q T < 2 GeV/c, because for larger cutoffs the asymmetry is diluted. 


2. A new high-energy parametrization 

As already mentioned above, this new parametrization is inspired to the one of Ref. [37]. There, it was assumed that 
the transverse momentum of the detected pion in the SIDIS process was entirely due to the transverse-momentum 
dependence in the Sivers function. No transverse momenta are contributed by other terms in the factorization 
formula. In this perspective, this approach can be considered as a limiting case of the approach of Ref. [38], based 
on Gaussian ansatze for the p T dependence both in the distribution and fragmentation functions (see Ref. [40] for a 
more detailed discussion). In Ref. [37], no further assumption was made but the p T distribution was integrated out. 
As a result, the SSA for the Sivers effect in SIDIS was expressed in terms of ” ^-moments” of the Sivers function, 
which were parametrized in terms of the w-quark distribution /“(x) and flavor-dependent normalizations S u ,Sd, to 
be determined by a fit to the new HERMES data [ 8 ]. Also in this case the normalizations turn out to have opposite 
sign, and the y 2 per degree of freedom is 1.2. The recent COMPASS data [10] were not included in the fit, but a 
direct comparison show a qualitative good agreement. 

Here, we retain the x dependence of this approach, but we introduce a different flavor-dependent normalization and 
an explicit p T dependence that are bound to the shape of the recent RHIC data on ppT —> wX at y/s = 200 GeV [4]. 
The expression adopted is 


fir (x,Pt) = N f x(l-x) 


= Nfx( 1 — x) 


M 2 pIpt 

(p 2 T +ir 

M 2 pIpt 


f{ {x,Pt) 

1 


(Pt + 1 1 ) 2 Vp 2 ) 

where po = 2 GeV/c. Following the same arguments of previous section, we get 


e ~pt/(pt) ff (a) 


/ 2 p 0 q T \ 2 E/ e/ N f h(xi;f /p) fi(x 2 ; / /p ] ) 
\Q 2 +Po) E/ e } h(xi-, f /p) f{ (x 2 ; f /p) + (1 <-> 2) 


(15) 


(16) 


1 In Eq. (14), the factor in front of the flavor-dependent term should read 1/18 because of the symmetry operation in the denominator 
of Eq. (13). However, as it is shown in the next Sec. IV A the SSA is not suppressed only in the (xi < 0, X 2 > 0) region of the phase 
space, which corresponds to take the dominant part of just the first term in the denominator of Eq. (13). 



Again, we can further simplify the expression introducing the flavor average product (/i(xi)) (fi(x 2 )}, which leads 
to 2 


C 4 


\S 2T \ X2 (1 - X 2 ) 


2 Po q T \ 2 8 N u + N d 
Qt + Po) 9 


(17) 


The q T shape is different from Eq. (14) and the peak position is shifted at much larger values. This is in agreement 
with a similar analysis of the azimuthal asymmetry of the unpolarized Drell-Yan data (the violation of the Lam- 
Tung sum rule [41]). But, more specifically, it is induced by the observed x F — q T correlation in the RHIC data for 
pp^ —> 7 tX, when it is assumed that the SSA is entirely due to the Sivers mechanism; this suggests that the maximum 
asymmetry is reached in the upper valence region such that x F ss x 2 ~ (<? r )/5 [4]. We have conveniently modified 
the cutoffs such that for this parametrization the sampled distribution is 1 < q T < 3 GeV/c. In this case, the peak 
asymmetry is reached for x 2 ~ 0.5 (see next Sec. IV). Contrary to the other Lorentzian and Gaussian distributions 
adopted so far, the p T distribution of Eq. (15) cannot be normalized in the usual way because of lack of convergence 
at large q T . Rather, it is built as to assume the value 1 on its peak position q T = po- Finally, the flavor dependence 
of the normalization is kept as simple as in Ref. [38], namely N u = —Nd = 0.7. The sign, positive for u quark and 
negative for the d, already takes into account the predicted sign change of fi T from Drell-Yan to SIDIS, where the 
opposite flavor dependence of the sign was obtained [36-38]. 


IV. RESULTS OF THE MONTE CARLO SIMULATIONS 

In this Section, we present results for Monte Carlo simulations of the Sivers effect in the Drell-Yan process pp'' —> 
p + X using input from the previous Sec. Ill 1 and III 2. The goal is to explore the sensitivity of the simulated 
asymmetry to different input parametrizations (i.e. to the input theoretical uncertainty), as well as to directly verify, 
within the reached statistical accuracy, the predicted sign change of the Sivers function between SIDIS and Drell- 
Yan [31]. The collision is considered at the cm energy y/s = 200 GeV, at which RHIC is presently taking data. We 
use a conservative dilution factor 0.5 for the proton polarization, even if a foreseen upgrade of the superconducting 
Siberian snake will allow RHIC to run in the future with stable 70% transverse polarization [55]. We select two 
different ranges for the lepton invariant mass: 4 < M < 9 GeV and 12 < M < 40 GeV. In this way, we avoid overlaps 
with the resonance regions of the cc and bb quarkonium systems. At the same time, the theoretical analysis based on 
the leading-twist cross section (2) should be well established, since higher-twist effects can be classified according to 
powers of M p /M, where M p is the proton mass. Moreover, at y/s = 200 GeV also the QCD corrections beyond tree 
level should be suppressed, which justifies the approximations described in Sec. II. In the Monte Carlo, the events 
are sorted according to the cross section (5), supplemented by Eqs. (6)-(8), while the asymmetry related to the Sivers 
effect is simulated by Eqs. (9), (10), (14) and (17). In particular, the events are divided in two groups, one for positive 
values (U) of sin((/>—</>s 2 ) in Eq. (10), and another one for negative values (D), and taking the ratio (U — D)/(U + D). 
Data are accumulated only in the x 2 bins of the polarized proton, i.e. they are summed over in the x\ bins for the 
unpolarized proton, in the transverse momentum q T of the muon pair and in their zenithal orientation 8. Proper cuts 
are applied to the q T distribution according to the different input parametrization of the Sivers function: for the case 
of Sec. Ill 1, 0.1 < q T < 2 GeV/c; for the case of Sec. Ill 2, 1 < q T < 3 GeV/c. In this way, the ratio between the 
absolute sizes of the asymmetry and the statistical errors is optimized for each choice. The resulting (q T ) is ~ 1.8 
GeV/c, in fair agreement with the one experimentally explored at RHIC [4], Contrary to Refs. [42, 46], there is no 
need to introduce cuts in the 8 distribution because of the (1 + cos 2 8) term in Eq. (10). We have considered two 
initial different samples of 20 000 and 100 000 events. Statistical errors for (U — D)/(U + D) are obtained by making 
10 independent repetitions of the simulation for each individual case, and then calculating for each x 2 bin the average 
asymmetry value and the variance. We checked that 10 repetitions are a reasonable threshold to have stable numbers, 
since the results do not change significantly when increasing the number of repetitions beyond 6. 


A. Properties of the phase space 

In Fig. 2, the scatter plot (in the fractional momenta x\ = x ', x 2 = x, of the annihilating partons) for 120 000 events 
of Drell-Yan muon pairs produced by proton-proton collisions at y/s = 200 GeV, is shown. The data are divided in 


2 An argument similar to the one in the previous footnote about Eq. (14), applies also here to Eq. (17). 
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FIG. 2: The scatter plot for 120 000 events of Drell-Yan muon pairs produced by proton-proton collisions at \fs = 200 GeV. 
The two bands in which data are grouped, correspond to two different ranges in the muon invariant mass: 4 < M < 9 GeV for 
the lower band, 12 < M < 40 GeV for the upper one. 


two different bands, corresponding to two different ranges in the muon pair invariant mass: 4 < M < 9 GeV for the 
lower band and 12 < M < 40 GeV for the upper one. In fact, hyperboles X\X 2 = const, are selected by fixed values 
of t = X\X 2 = M 2 / s. Since the elementary annihilation is assumed to proceed through a virtual photon, the cross 
section contains a term 1/M 2 ~ 1/r which populates the phase space at low values, while the upper right corner of 
Fig. 2 for r —» 1 is basically empty (also because the parton distributions vanish for X \/2 —► 1). Therefore, within 
each band events tend to accumulate to the lowest possible M in the considered range, which means that they try to 
align along the hyperbole with lowest possible values of X\ and X 2 ■ Moreover, for the same reason the lower band is 
much more dense than the other one: 95% of the events correspond to the 4 < M < 9 GeV range. This is why we 
consider also this case, because the much higher statistics can be traded for the questionable neglect of higher-twist 
contributions with respect to the higher 12 < M < 40 GeV range. 

The scatter plot of Fig. 2 contains two main differences with respect to our previous analysis of Ref. [42] for the GSI 
setup, where the pj/ —> process was considered at ^fs = 14 GeV. First of all, here the cm energy is higher 

by one order of magnitude, which means that Xi / 2 values lower by one order of magnitude are explored, typically 
( 24 / 2 ) ~ 0.01. Secondly, this is emphasized by the fact that in pp collisions at least one of the two annihilating partons 
comes from the proton sea distribution, which is peaked at very low x values. The importance of parton momenta 
below the valence region is potentially relevant to the theoretical models. We already mentioned in Sec. Ill 1 that 
the parametrization of Eq. (12) leads to a persistence of the Sivers effect in this range, contrary to the other one of 
Ecp (15). Anyway, both choices lead to a vanishing effect for X 2 —> 0. The dominance of the low X 2 portion of phase 
space is evident in the histograms for the event distributions displayed in the next Sec. IV, where the first bin [0,0.1] 
contains more than 50% of events on average. We can also estimate the expected position of the peak density, assuming 
that the bidimensional event distribution N(x\,X 2 ) is dominated by the 1 /r factor associated to the elementary qq 
fusion into a virtual photon. In this case, N(x 1 , 2 : 2 ) = l/{x\X 2 ) for t > M 2 %in /s, and 0 otherwise. Consequently, the 
X\ -integrated distribution has the form \og(x2s/M/ nin )/x2 and reaches its peak value for log(x 2 s/M^ lin ) = 1, i.e. for 
X 2 = eM'/ nin /s rts 2,M/ nin /s ~ 0.01 for yfs = 200 GeV and M min = 12 GeV, in agreement with the average values 
produced by the Monte Carlo. Elsewhere, the integrated distribution behaves like 1 / 2 : 2 , as it can visually be checked 
in the histograms of next Sec. IV. 

The 1/r mechanism in the random generation can artificially suppress the asymmetry irrespective of the size of 
the Sivers function itself. In fact, let us consider the 12 < M < 40 GeV range, which has a lower event rate; we 
distinguish four different slices of phase space: 

• xi,X 2 > 0.1 : this part covers 99% of the phase space, but it contains only 0.5% of the total number of events; 
it corresponds to higher M values (M > 20 GeV) and it is suppressed by the 1/M 2 mechanism; moreover, in 
F(x 1, X 2 ) of Eq. (8) the annihilating antiquark with flavor / is picked up from the sea distribution of one of the 
two protons at large x; 

• X\,X 2 < 0.1 : this part covers < 1 % of the phase space, but it contains 20 % of events; in fact, it corresponds to 
lower M values (emphasized by the 1/M 2 mechanism) and F[x 1,2:2) is dominated by the sea distributions in 



10 


TABLE II: Total cross sections for Drell-Yan pp collisions at y/s = 200 GeV and for various invariant masses of the muon pair, 
at the given luminosity 10 32 cm -2 sec -1 . 


fi{x) from 

M (GeV) 

a pp (nb) 

rates (events/month) 

Ref. [14] 

44-9 

1.2 

2.5 x 10 5 

Ref. [14] 

12 -4 40 

0.1 

2.5 x 10 4 

Ref. [56] 

4^9 

7 

1.5 x 10 s 

Ref. [56] 

12 4-40 

0.4 

10 5 


both protons, that are enhanced at small x: however, the SSA is suppressed because fyp/fi —> 0 for X 2 —> 0; 

• x\ > 0.1, X 2 < 0.1 : this part covers again < 1% of the phase space, but it contains 40% of the events; it 
is less favoured by the 1/M 2 mechanism with respect to the previous case, but F(x 1 , 2 : 2 ) contains the term 

f /p) fi(x 2 ', f /p), that is dominant in this slice of phase space;_however, for the very same reason the SSA 
is suppressed because it is approximately driven by f/p{x 2 )/fi(x 2 ', f /p) —> 0 for X 2 —> 0; 

• X\ < 0.1, X 2 > 0.1 : all previous arguments apply also here upon the X\ <-> X 2 exchange but for the SSA, which 
is driven by fi T (x 2 )/fi(x 2 ]f/p) and, therefore, it is not suppressed for X 2 > 0.1. 

In summary, irrespective of the size of the Sivers function, the magnitude of the corresponding SSA is suppressed in 
all parts of phase space but in the region x\ < 0.1, X 2 > 0.1, dominated by the sea partons of the unpolarized proton 
and by the valence partons of the polarized one. Elsewhere, the 1/M 2 mechanism induces a dominance of the sea 
partons, that acts as an effective dilution factor leading to a waste of ~ 50% of the total number of events. 

B. Total cross section and event rates 

Using Eq. (8) with the parametrization for the parton distributions from Ref. [14], from our Monte Carlo we deduce 
a total cross section a pp = 0.1 nb for the pp —> p + X process at \/~s = 200 GeV and with invariant masses in the 
range 12 < M < 40 GeV, while we get a pp = 1.2 nb for the lower 4 < M < 9 GeV range. The results are quite 
sensitive to the parametrization of the parton distributions. We have recalculated the total cross sections with the 
more recent NNLO analysis of Ref. [56] and we get 0.4 nb and 7 nb, respectively. When changing parametrizations 
and, consequently, normalizations, we had to readjust the K factor accordingly; in order to reproduce the measured 
cross section at t/s « 16 GeV [14], we had to reduce it by 50%. It means that now I\ « 1 and the QCD corrections 
are mostly contained in the NNLO parametrization of the parton distributions. The sensitivity (and the related 
uncertainty) of this analysis to the input are sizeable, but do not alter the order of magnitude of the result. Since 
our goal is to estimate event rates by multiplying a pp with a given luminosity, we are confident that the results are 
realistic and reliable. For RHIC, a luminosity of 10 32 cm -2 sec -1 or higher is foreseen [48]. This means, for example, 
that at least 250 000 Drell-Yan events/month (and up to 7 times more) could be collected with this luminosity and 
muon pair invariant masses in the 4 < M < 9 GeV range. A list of the combinations here explored is given in Tab. II 
(for a more comprehensive analysis see Ref. [57]). 

It is also interesting to compare with the antiproton-proton Drell-Yan collision. In general, we expect that the lower 
the ( x ), the more the cross sections are dominated by sea parton distributions, the closer the ratio <J pp /a pp approaches 
1. By updating our previous results [42] at the present energy and viceversa, we get 

^(v/s = 200; 4 < M < 9) = 2 = 200; 12 < M < 40) = 4 ^-{^s = 14; 4 < M < 9) = 40 . (18) 

Cpp &pp ®pp 

There are two ways to lower the range of x or, equivalently, r: decreasing the invariant mass M, or increasing the cm 
energy -^/s. Indeed, the results show that for this trend the ratio approaches 1. When increasing ^fs at a given M 
range, for example, the depletion of the ratio implies also that a pp increases. This curious result of an increasing cross 
section with energy can be explained by recalling that a shift to smaller Xi,X 2 , makes F( £ 1 , 2 : 2 ) in Eq. (8) dominated 
by the sea parton distributions, which are large at very small parton momenta. 

C. Single-spin asymmetries 

In Fig. 3, the sample of 20 000 Drell-Yan events for the pp t —> p + X reaction at yfs = 200 GeV is displayed for 
muon invariant mass in the 12 < M < 40 GeV range. Results are reported in X 2 bins excluding the upper boundary 
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FIG. 3: The sample of 20 000 Drcll-Yan events for the pp' —> pTX reaction at \fs = 200 GeV, 12 < M < 40 GeV, and 
1 < qr < 3 GeV/c, using the parametrization of Eq. (15) (see text). Upper panel: for each bin in the parton momentum X 2 
inside p\ the darker histogram collects events with positive sin (cf) — (j>s 2 ) (G), the superimposed lighter histogram collects the 
negative ones ( D ). Lower panel: the asymmetry (U — D)/(U + D); upward triangles for N u > 0 in Eq. (15) corresponds to 
a sign change in the Sivers function from SIDIS to Drell-Yan processes; squares for N u < 0. Statistical error bars from 10 
independent repetitions of the simulation. Continuous lines are drawn to guide the eye. 


X 2 > 0.8, which is scarcely or not at all populated, according to the upper band in Fig. 2. Events are accumulated 
according to Eq. (17) based on the parametrization (15) of the Sivers function; consequently, the transverse momentum 
distribution is constrained by 1 < q T < 3 GeV/c. In the upper panel, for each bin two groups of events are stored, 
one corresponding to positive values of sin(</> — (j>s 2 ) in Eq. (10) (represented by the darker histogram), and one for 
negative values (superimposed lighter histogram). In the lower panel, the asymmetry (U — D)/(U + D) is shown 
between the positive ( U ) and negative ( D) values. Average values of the asymmetry and (statistical) error bars 
are obtained by 10 independent repetitions of the simulation. The upward triangles indicate the results assuming a 
positive normalization for the quark u in Eq. (15), which already takes into account the predicted sign change of fi T 
from Drell-Yan to SIDIS [31] with respect to recent parametrizations of SIDIS data [36-38]. For sake of comparison, 
the squares illustrate the opposite results that one would obtain ignoring such prediction. 

From the upper panel of Fig. 3, we deduce that the assumed elementary qq —> 7 * mechanism indeed populates the 
phase space for the lowest possible r, with more than 50% of the events in the 0 < X 2 < 0.1 bin, leaving a ~ I/X 2 
distribution outside. In the lower panel, correspondingly, the error bars are small for X 2 < 0.5 and allow for a clean 
reconstruction of the asymmetry shape and, more importantly, for a conclusive test of the predicted sign change in 
fyp. With the considered sample of 20 000 events, the same conclusion is not possible using the parametrization of 
Eq. (12), because the asymmetry produced by Eq. (14) is too small. More quantitatively, half of the displayed error 
bar representes the variance AA(x- 2 ) for the asymmetry A = (U — D)/(U + D) in each X 2 bin. The results in the lower 
panel of Fig. 3 can be approximated by the relation A A « 0.05x2. The asymmetry is statistically not compatible 
with zero if A(x 2 ) > A A = 0.05x2 for the considered X 2 range. With 20 000 events, this condition is fulfilled only 
by the parametrization (15), but not by the one in Eq. (12). Finally, from Tab. II we deduce that an hypothetical 
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experiment in these kinematic conditions should run from one week to, at most, almost one month in order to reach 
the indicated statistical error bars. 


30000 
25000 
20000 
15000 
10000 
5000 
0 

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 

0.08 
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- 0.06 
- 0.08 
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X 2 

FIG. 4: Same as in Fig. 3 but considering 100 000 events using the parametrization of Eq. (12) and 0.1 < qr < 2 GeV/c (see 
text). Lower panel: upward triangles for N u > 0 in Eq. (14), squares for N u < 0. 




events 



x 2 x 2 

FIG. 5: The sample of 100 000 Drcll-Yan events for the same reaction in the same conditions and notations as in Fig. 3 but for 
4 < M < 9 GeV. Left panel: parametrization of Eq. (15); right panel: parametrization of Eq. (12) (see text). 

In Fig. 4, we consider the same kinematic conditions of the previous Fig. 3 but with a sample of 100 000 events, 
which would require at least one month of running time, or four months in the most disfavoured conditions. We 
employ the parametrization of Eq. (12). We recall that the q T distribution is now integrated in the range 0.1 < q T < 2 
GeV/c with a resulting lower (q%) (see previous Sec. Ill 1). Notations in the figure are the same as in the previous 
Fig. 3. In particular, in the lower panel upward triangles identify the results from Eq. (14) with N u > 0 from Tab. I; 
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squared points refer to the opposite choice. Therefore, we conclude that the parametrization (12) demands for a 
much higher statistics in order to get a clear nonvanishing shape, because it produces an overall smaller asymmetry. 
Still, a definite answer is possible about the sign change prediction of fyp if the statistical sample meets the required 
conditions. 

In Fig. 5, we show just the histogram of collected 100 000 events in the same conditions and notations as before but 
for the lower 4 < M < 9 GeV range. From Tab. II, we note that the necessary running time is 24-3 times shorter than 
the one for collecting 20 000 events at the higher 12 < M < 40 GeV range, depending on the parametrization chosen 
for f\. Again, this is due to the already mentioned 1/r factor of the elementary qq —> 7 * mechanism, which privileges 
lower r. For the very same reason, an even larger portion of events (77%) is contained in the first 0 < X 2 <0.1 bin, 
while the remaining 23% is distributed for x 2 > 0.1 approximately as l/x 2 - 


asymmetry 



FIG. 6: The asymmetry (U — D)/(U + D) corresponding to the histograms of Fig. 5, where U identifies the darker histograms 
and D the superimposed lighter ones (see text). Upward triangles for the parametrization of Eq. (15) with N u > 0; squares for 
N u < 0. Open upward triangles for the parametrization of Eq. (12) with N u > 0; open squares for N u < 0. 

In Fig. 6, we plot the spin asymmetries corresponding to the histograms of the previous Fig. 5. Notations are the 
following. Upward triangles correspond to the histograms in the left panel of Fig. 5, i.e. to the parametrization of 
Eq. (15); the flavor dependent normalization is N u > 0, according to the properties of f^ T as it is extracted from 
SIDIS data. The corresponding choice N u < 0 is represented by the squares. The open upward triangles and squares 
show the results for the other parametrization Eq. (12). The accumulation of events for very low X 2 values, which 
is evident in Fig. 5, here reflects in very tiny error bars for the same bins, allowing to clearly distinguish the two 
different parametrizations for 0 < X 2 < 0.5. In this range, a measurement of the Sivers effect for invariant masses as 
low as 4 < M < 9 GeV, allows to reduce the theoretical uncertainties in very few days of running time (from 2 to 12, 
depending on the parametrization of the unpolarized parton distributions). Moreover, and most important, in the 
same X 2 range the statistical accuracy is sufficient to directly test the predicted sign change of fyj* [31], irrespective 
of the uncertainty in the theoretical input. 

For 0.5 < X 2 , the width of the error bars does not always allow for such analysis, since the variance grows approx¬ 
imately as A A ss 0.04x2- We already noted that for 12 < M < 40 GeV approximately half of the 20 000 events lie 
in the valence range 0.1 < X 2 < 0.7 (see the upper panel in Fig. 3), while for 4 < M < 9 GeV 23 000 events out 
of 100 000 are found in the same range (see Fig. 5). If we assume that the size of the error bars is proportional to 
1/VN, with N the number of events, then, in the valence region for a given parametrization, the size of the error 
bars for the higher M range should approximately scale as \J 23/10 ss 1.5 with respect to the one for the lower M 
range. It is easy to check from the upward triangles in Fig. 3 and Fig. 6 that for the parametrization (15) indeed this 
approximate relation is verified. But from Tab. II, we deduce also that in the same running time necessary to collect 
20 000 events at 12 < M < 40 GeV it is possible to collect 200 000 4- 300 000 events at 4 < M < 9 GeV, depending on 
the chosen parametrization. In the valence range, this means 46 000 4- 70 000 events (the 23% of the total, as before), 
which induces a reduction factor 2 4- 2.5 in the size of the error bars. Therefore, we can put a sort of ’’normalization” 
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for the approximate behaviour of the size of the variance, by guessing that 

A . nnr 1 20 000 
A A(x 2 ) « 0.05 x 2 y N 

\ t / . nnt 1 100 000 

A A(x 2 ) ~ 0.04x2 Y —— 

For the lowest 0 < x 2 < 0.1 bin the coefficient is 0.008 and 0.004, respectively, and the scale factor in the size of the 
error bars gets amplified. 


12 < M < 40 GeV 

0.1 < x 2 < 0.7 

(19) 

4 < M < 9 GeV 

0.1 < x 2 < 0.7 . 

(20) 




0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 

X 2 


FIG. 7: The asymmetry corresponding to the open upward triangles and open squares in Fig. 6, is shown in the 0 < X 2 < 0.2 
range with a finer binning. 

In Fig. 7, we show a finer binning of the range 0 < x 2 <0.2 for the asymmetry obtained with the parametrization (12) 
for 100 000 events in the 4 < M < 9 GeV case. It is a closer view of the very tiny error bars of open upward triangles 
and squares in Fig. 6. This is a peculiar feature of this parametrization, which emphasizes the role of very low x 2 
through the parameters in Tab. I, as it is discussed in Sec. Ill 1. 


V. CONCLUSIONS 

In this paper, we have concentrated on the investigation of the spin structure of the proton using the single- 
polarized Drell-Yan process pp^ —► p, + p,~X. At leading twist, the cross section contains several terms that lead to 
an asymmetric distribution of the final muon pair in its azimuthal angle (f> with respect to the production plane. In 
previous papers [42, 46], we considered those terms involving the transversity distribution hi and the Boer-Mulders 
function hj~ [41], which is believed to be responsible for the very well known violation of the Lam-Tung sum rule in the 
unpolarized Drell-Yan data [12-14], There, we set up a Monte Carlo to numerically simulate (polarized) antiproton- 
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proton Drell-Yan collisions to study the best kinematic conditions for the HESR at GSI [43, 44] that allow to extract 
unambiguous information on the target parton distributions. 

Here, we have followed the same approach to isolate the contribution of the term involving the Sivers function /'/y, 
a ’’naive T-odd” TMD partonic density that describes how the distribution of unpolarized quarks is distorted by the 
transverse polarization of the parent hadron. As such, //y contains unsuppressed information on the orbital motion 
of hidden confined partons; better, it contains information on their spatial distribution [29], and it offers a natural link 
between microscopic properties of confined elementary constituents and hadronic measurable quantities, such as the 
nucleon anomalous magnetic moment [30]. Factorization theorems for TMD distribution and fragmentation functions 
indicate that //y is universal modulo a sign change (when switching from the SIDIS to the Drell-Yan process), due 
to a peculiar feature under the time-reversal operation of the gauge link operator required to make its definition 
color-gauge invariant [31]. 

Recently, very precise data for SSA involving /fly (the Sivers effect) have been obtained for the SIDIS process on 
transversely polarized protons [ 8 ]. This allowed for more realistic parametrizations of fi T [36-38], that have been 
used then to make predictions for SSA in proton-proton collisions at RHIC (for a comparison among the various 
approaches, see also Ref. [40]). 

Here, we have numerically simulated the Sivers effect for the pp^ —> p + p~X process at y/s = 200 GeV including the 
foreseen upgrade in the RHIC luminosity (RHIC II). The goal is to explore the sensitivity of the simulated asymmetry 
to different input parametrizations, as well as to directly verify, within the reached statistical accuracy, the predicted 
sign change in the universality properties of the Sivers function. Therefore, we have employed the parametrization 
of Ref. [36] and a new high-energy parametrization of //y, whose flavor-dependent normalization and p T distribution 
are constrained by recent RHIC data on SSA for the pp^ —> ttX process at yfs = 200 GeV [4]. The main difference is 
that the former, fitted to the SIDIS data of Ref. [ 8 ], displays an emphasized relative importance of the unfavoured d 
quark, and it gives an average transverse momentum (q T ) of the lepton pair much lower than the latter. Consistently, 
we have built spin asymmetries by integrating the q T distribution with adequate cutoffs, namely 0.1 < q T < 2 GeV/c 
for the former parametrization, and 1 < q T < 3 GeV/c for the latter. Results have been presented as binned in the 
parton momenta X 2 of the polarized proton, i.e. by integrating also upon the antiparton partner momenta Xi and the 
zenithal muon pair distribution 9 with no further cuts. 

Sorted events are divided in two groups, corresponding to opposite azimuthal orientations of the muon pair with 
respect to the reaction plane (conventionally indicated with U and D ), and the asymmetry (U — D)/(U + D) has 
been considered. Two different samples of 20 000 and 100 000 events have been selected, and statistical errors for 
(U — D)/(U + D) have been obtained by making 10 independent repetitions of the simulation for each individual case 
and, then, calculating for each X 2 bin the average asymmetry and the variance. We have considered two different 
ranges of muon pair invariant mass, namely 4 < M < 9 GeV and 12 < M < 40 GeV. In this way, we avoid overlaps 
with the resonance regions of the cc and bb quarkonium systems, and we can safely assume that the elementary 
annihilation proceeds through the qq —> 7 * mechanism. In particular, the Monte Carlo is based on the corresponding 
leading-twist cross section. In the higher mass range, this theoretical analysis appears well established, since higher 
twists may be suppressed as M p /M, where M p is the proton mass. More questionable is the case of the lower M 
range, but this uncertainty can be traded for the much higher statistics that can be reached because of the 1/M 2 
contribution of the 7 * propagator. Indeed, approximately 95% of the events fall in the 4 < M < 9 GeV range with a 
significant reduction of the running time necessary to reach a predefined statistical accuracy. 

More generally, a very small portion of the phase space, corresponding to the lowest possible values of r = 2727 = 
M 2 /s, contains most of the events, also because of the high cm square energy s; typically, ( 27 / 2 ) ~ 0.01. This is 
emphasized by the fact that in pp collisions at least one of the two annihilating partons comes from the proton sea 
distribution, which is peaked at very low x values. The importance of parton momenta below the valence region is 
potentially relevant to the theoretical models. In our case, it turns out that the parametrization of Ref. [36] gives 
small asymmetries because the emphasized unfavoured dd annihilation is statistically suppressed. As a consequence, 
in the 12 < M < 40 GeV range 20 000 events are not sufficient to produce a Sivers effect that is not statistically 
consistent with zero. On the contrary, with the high-energy parametrization described in this paper this analysis is 
possible in the range 0 < 27 < 0.5, and a direct test of the ’’universal” sign change of fy T can be unambiguously 
performed. 

A more favourable situation is encountered with lower M values, say in the 4 < M < 9 GeV range. The much 
higher statistics significantly reduces the error bars and allows for a very short running time, typically as short as 
few days to collect the 100 000 events simulated here at the foreseen luminosity of 10 32 cm -2 sec -1 , irrespectively of 
the choice of input parametrizations for f^p. Moreover, we observe that in the 0 < 27 < 0.5 range the two choices 
give two clearly distinct asymmetries, and for each case the results with opposite signs in the Sivers function can be 
unambiguously separated. Remarkably, we stress that this analysis holds also at very low 27 values, typically as low 
as 0.01, which are relevant at RHIC. 

In conclusion, at the foreseen luminosity of 10 32 cm -2 sec -1 (RHIC II) and with a dilution factor 0.5, in few days 
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RHIC can collect 100 000 Drell-Yan events for the process pp^ —> p + X at ^fs = 200 GeV and with muon pair 
invariant masses in the 4 < M < 9 GeV range. From the measured Sivers effect, it should be possible to extract 
information on the x structure of the Sivers function in the 0 < x < 0.5 range, particularly also at very low x ~ 0.01, 
as well as to test its ’’universal” sign change predicted in Ref. [31]. At higher M values, like 12 < M < 40 GeV, 
the situation is theoretically more favourable because the higher twists are suppressed as 1/M. However, the lower 
density in the phase space makes the running time much longer. It is still possible to perform the previous analysis 
with 100 000 events and to come to the same conclusions, but at the price of taking data for some months. With a 
reduced sample, e.g. of 20 000 events, this time is shortened to few weeks, but our analysis was possible for only one 
of the chosen parametrizations, the other one giving results compatible with zero. 
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